home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 39 / Amiga Format CD39 (1999-04-13)(Future Publishing)(GB)[!][issue 1999-05].iso / -seriously_amiga- / graphics / mpeg2decode / src / mpegaudio / subs.c < prev   
C/C++ Source or Header  |  1999-03-02  |  7KB  |  167 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. subs.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@3d.enet.dec.com             *
  10.  *                                                                    *
  11.  * VERSION 3.9                                                        *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  15.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  16.  * 7/10/91  Earle Jennings      Ported to MsDos from Macintosh        *
  17.  *                              Replacement of one float with FLOAT   *
  18.  * 2/11/92  W. Joseph Carter    Added type casting to memset() args.  *
  19.  * 4/27/92  Masahiro Iwadare    Added 256 point version for Layer III *
  20.  **********************************************************************/
  21.  
  22. #include "common.h"
  23. #include "encoder.h"
  24.  
  25. /*****************************************************************************
  26.  ************************** Start of Subroutines *****************************
  27.  *****************************************************************************/
  28.  
  29. /*****************************************************************************
  30.  * FFT computes fast fourier transform of BLKSIZE samples of data            *
  31.  *   uses decimation-in-frequency algorithm described in "Digital            *
  32.  *   Signal Processing" by Oppenheim and Schafer, refer to pages 304         *
  33.  *   (flow graph) and 330-332 (Fortran program in problem 5)                 *
  34.  *   to get the inverse fft, change line 20 from                             *
  35.  *                 w_imag[L] = -sin(PI/le1);                                 *
  36.  *                          to                                               *
  37.  *                 w_imag[L] = sin(PI/le1);                                  *
  38.  *                                                                           *
  39.  *   required constants:                                                     *
  40.  *         #define      PI          3.14159265358979                         *
  41.  *         #define      BLKSIZE     1024                                     *
  42.  *         #define      LOGBLKSIZE  10                                       *
  43.  *         #define      BLKSIZE_S   256                                      *
  44.  *         #define      LOGBLKSIZE_S 8                                       *
  45.  *                                                                           *
  46.  *****************************************************************************/
  47. #define      BLKSIZE_S   256
  48. #define      LOGBLKSIZE_S 8
  49.  
  50. #ifdef STORM
  51. void fft(FLOAT x_real[BLKSIZE],FLOAT x_imag[BLKSIZE], FLOAT energy[BLKSIZE], FLOAT phi[BLKSIZE], int N)
  52. #else
  53. void fft(x_real,x_imag, energy, phi, N)
  54. FLOAT x_real[BLKSIZE], x_imag[BLKSIZE], energy[BLKSIZE], phi[BLKSIZE];
  55. int N;
  56. #endif
  57. {
  58.  int     M,MM1;
  59.  static int     init=0;
  60.  int     NV2, NM1, MP;
  61.  static double  w_real[2][LOGBLKSIZE], w_imag[2][LOGBLKSIZE];
  62.  int            i,j,k,L;
  63.  int            ip, le,le1;
  64.  double         t_real, t_imag, u_real, u_imag;
  65.  
  66.  if(init==0) {
  67.     memset((char *) w_real, 0, sizeof(w_real));  /* preset statics to 0 */
  68.     memset((char *) w_imag, 0, sizeof(w_imag));  /* preset statics to 0 */
  69.     M = LOGBLKSIZE;
  70.     for(L=0; L<M; L++){
  71.        le = 1 << (M-L);
  72.        le1 = le >> 1;
  73.        w_real[0][L] = cos(PI/le1);
  74.        w_imag[0][L] = -sin(PI/le1);
  75.     }          
  76.     M = LOGBLKSIZE_S;
  77.     for(L=0; L<M; L++){
  78.        le = 1 << (M-L);
  79.        le1 = le >> 1;
  80.        w_real[1][L] = cos(PI/le1);
  81.        w_imag[1][L] = -sin(PI/le1);
  82.     }          
  83.     init++;
  84.  }
  85.  switch(N) {
  86.     case BLKSIZE:
  87.             M = LOGBLKSIZE;
  88.             MP = 0;
  89.             break;
  90.     case BLKSIZE_S:
  91.             M = LOGBLKSIZE_S;
  92.             MP = 1;
  93.             break;
  94.     default:    printf("Error: Bad FFT Size in subs.c\n");
  95.             exit(-1);
  96.  }
  97.  MM1 = M-1;
  98.  NV2 = N >> 1;
  99.  NM1 = N - 1;
  100.  for(L=0; L<MM1; L++){
  101.     le = 1 << (M-L);
  102.     le1 = le >> 1;
  103.     u_real = 1;
  104.     u_imag = 0;
  105.     for(j=0; j<le1; j++){
  106.        for(i=j; i<N; i+=le){
  107.           ip = i + le1;
  108.           t_real = x_real[i] + x_real[ip];
  109.           t_imag = x_imag[i] + x_imag[ip];
  110.           x_real[ip] = x_real[i] - x_real[ip];
  111.           x_imag[ip] = x_imag[i] - x_imag[ip];
  112.           x_real[i] = t_real;
  113.           x_imag[i] = t_imag;
  114.           t_real = x_real[ip];
  115.           x_real[ip] = x_real[ip]*u_real - x_imag[ip]*u_imag;
  116.           x_imag[ip] = x_imag[ip]*u_real + t_real*u_imag;
  117.        }
  118.        t_real = u_real;
  119.        u_real = u_real*w_real[MP][L] - u_imag*w_imag[MP][L];
  120.        u_imag = u_imag*w_real[MP][L] + t_real*w_imag[MP][L];
  121.     }
  122.  }
  123.  /* special case: L = M-1; all Wn = 1 */
  124.  for(i=0; i<N; i+=2){
  125.     ip = i + 1;
  126.     t_real = x_real[i] + x_real[ip];
  127.     t_imag = x_imag[i] + x_imag[ip];
  128.     x_real[ip] = x_real[i] - x_real[ip];
  129.     x_imag[ip] = x_imag[i] - x_imag[ip];
  130.     x_real[i] = t_real;
  131.     x_imag[i] = t_imag;
  132.     energy[i] = x_real[i]*x_real[i] + x_imag[i]*x_imag[i];
  133.     if(energy[i] <= 0.0005){phi[i] = 0;energy[i] = 0.0005;}
  134.     else phi[i] = atan2((double) x_imag[i],(double) x_real[i]);
  135.     energy[ip] = x_real[ip]*x_real[ip] + x_imag[ip]*x_imag[ip];
  136.     if(energy[ip] == 0)phi[ip] = 0;
  137.     else phi[ip] = atan2((double) x_imag[ip],(double) x_real[ip]);
  138.  }
  139.  /* this section reorders the data to the correct ordering */
  140.  j = 0;
  141.  for(i=0; i<NM1; i++){
  142.     if(i<j){
  143. /* use this section only if you need the FFT in complex number form *
  144.  * (and in the correct ordering)                                    */
  145.        t_real = x_real[j];
  146.        t_imag = x_imag[j];
  147.        x_real[j] = x_real[i];
  148.        x_imag[j] = x_imag[i];
  149.        x_real[i] = t_real;
  150.        x_imag[i] = t_imag;
  151. /* reorder the energy and phase, phi                                        */
  152.        t_real = energy[j];
  153.        energy[j] = energy[i];
  154.        energy[i] = t_real;
  155.        t_real = phi[j];
  156.        phi[j] = phi[i];
  157.        phi[i] = t_real;
  158.     }
  159.     k=NV2;
  160.     while(k<=j){
  161.        j = j-k;
  162.        k = k >> 1;
  163.     }
  164.     j = j+k;
  165.  }
  166. }
  167.